/*
 *	Trabalho 2 - Hospital em ICMCTown
 *	Programacao Concorrente - SCC0143 - 2014
 *
 *	Henrique Teixeira Arroyo	6419241
 *	Henry D. Leão				5881632
 *	Júlio Zanzini Polzato		7288971
 * 
 *  - Versão paralela do algoritmo utilizado
 *
 *  A paralelizacao foi realizada utilizando MPI e openMP
 * 
 */
 
#include <omp.h>
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h> 
#include <string.h>

#define NUM_THREADS 1		// Numero de threads a serem executadas em cada nó do cluster
#define INF 1000000			// Valor que representa nenhuma conexao entre dois vertices

float L;					//Raio do centro absoluto
int **D;					// Matrizes a serem utilizadas
int localDone;				// Valor que indica o fim do calculo da matriz de distancias

int	V,						// numero de nos
	E,						// numero de arestas
	myrank, 				// rank do processo
	P, 						// numero de processadores
	ns;						// numero de vertices no conjunto S

int	**edge,					// guarda as arestas do grafo edge[i][0]: nó fonte da aresta i, edge[i][1]: nó alvo da aresta i, edge[i][2]: peso da aresta i.
	*aw,					// posicao i em aw é a coordenada 'a' do par (a,b) do vertice w(i).
	*bw,					// posicao i em bw é a coordenada 'b' do par (a,b) do vertice w(i).
	*S;						// guarda valores dos vertices que permanecem apos a reducao.

int from,to;				// intervalo de iteracao de cada processo

//Funcao imprime a matriz 'M' de 'l' linhas e 'c' colunas
void print(int **M, int l, int c){
	int i,j;

	for(i=1;i<l;i++){
		for(j=0;j<c;j++)
			printf("%d ",M[i][j]);
		printf("\n");
	}
}

// Funcao ordena o vetor 'a' e realiza as mesmas operacoes no vetor 'b'
// os dois vetores devem ter dimensao N
void sort(int *a, int*b, int N){
	int xa, xb, i, j, aux;
	for(i=0;i<N;i++){
		xa = a[i];
		xb = b[i];
		j = i;
		while(j>0 && a[j-1] < xa){
			a[j] = a[j-1];
			b[j] = b[j-1];
			j--;
		}
		a[j] = xa;
		b[j] = xb;
	}
}

// Funcao executa a versao paralela do algoritmo de Floyd Warshal usando openMP
void runFloydWarshall(){
	int i, j, k;

	// Numero de linhas da matriz D enviadas para cada processo
	int task_rows = V/P;					
	// Numero de linhas a mais que sao enviadas para o ultimo processador quando V nao é divisivel por P
	int extra_rows = V%P;					

	// Calcula intervalo de iteracao de cada processo
	from = myrank*(task_rows);			
	to = from+task_rows;

	// Soma as linhas extras para o ultimo processo
	if(myrank==P-1)
		to+=extra_rows;

	int *displs = (int*)malloc(sizeof(int)*P);
	int *sendcount = (int*)malloc(sizeof(int)*P);
	// Calcula memoria a ser recolhida de cada processo
	for(i=0;i<P;i++){
		displs[i] = i*(task_rows)*V;		// offset a partir do qual a leitura e feita.
		sendcount[i] = (task_rows)*V;		// inteiros a serem lidos a partir do offset dado.
	}
	sendcount[P-1] += (extra_rows)*V;		// valores extras do ultimo processo.

	localDone = 0;
	// Inicia calculo da matriz de distancias
	while(!localDone){
		
		localDone = 1;
		// O loop é repetido até nenhum processo encontrar caminhos diferentes.
		for (i=from; i<to; i++){
			#pragma omp parallel for shared(D, V, i) private(j, k) reduction(&& : localDone)
			for (j=0; j<V; j++){
				for (k=0; k<V; k++){
					if(D[i][j] > D[i][k]+D[k][j]){
						D[i][j] = D[i][k]+D[k][j];
						localDone = 0;
					}
				}
			}
		}

		// Recolhe e redistribui os valores calculados
		MPI_Allreduce(&localDone, &localDone, 1, MPI_INT, MPI_BAND, MPI_COMM_WORLD);
		MPI_Allgatherv(D[from], (to-from)*V, MPI_INT, D[0], sendcount, displs, MPI_INT, MPI_COMM_WORLD);	
	}

	free(sendcount);
	free(displs);
}

// Funcao encontra o centro absoluto do grafo e retorna o seu raio na variavel global L.
void findAbsoluteCenter(){
	int i, j, b;
	float aL;

	// Numero de linhas da matriz D enviadas para cada processo
	int task_edges = E/P;					
	// Numero de linhas a mais que sao enviadas para o ultimo processador quando V nao é divisivel por P
	int extra_edges = E%P;					
	
	// Calcula intervalo de iteracao de cada processo
	from = myrank*(task_edges);			
	to = from+task_edges;
	
	// Soma as linhas extras para o ultimo processo
	if(myrank==P-1)
		to+=extra_edges;
	
	// Aloca variaveis
	aw = (int*)malloc(sizeof(int)*V);
	bw = (int*)malloc(sizeof(int)*V);	
 	S =  (int*)malloc(sizeof(int)*V);
	S[0] = 0;

	// Inicia algoritmo
	float l = INF;
	for(i=from;i<to;i++){
		// Encontra o par (a,b) para cada vertice
		#pragma omp parallel for shared (V,aw,bw,D) private(j)
		for(j=0;j<V;j++){
			aw[j] = D[edge[i][0]][j];
			bw[j] = D[edge[i][1]][j];
		}
		// Ordena conjunto de pontos
		sort(aw,bw,V);
		
		// Encontra vertices que entram no conjunto S apos a reducao.
		b = bw[0];
		ns = 1;
		for(j=1;j<V;j++){
			if(bw[j]>b){
				S[ns] = j;
				ns++;
				b = bw[j];
			}
		}
	
		// Calcula interseccoes (L) e guarda o menor valor que será o raio desejado.
		for(j=0;j<ns-1;j++){
			aL = ( edge[i][2] + bw[S[j]] + aw[S[j+1]] )*0.5;
			if(aL < l)
				l = aL;
		}		
	}
	
	// Reducao encontra o menor valor para o raio encontrado pelos processos criados.
	MPI_Reduce( &l, &L, 1, MPI_INT, MPI_MIN, 0, MPI_COMM_WORLD);
	
	free(aw);
	free(bw);
	free(S);
}

/*
 * Funcao aloca uma matriz de 'l' linhas e 'c' colunas de maneira contigua na memória.
 * (valores armazenados sequencialmente).
 */
int** allocMatrix(int l, int c){
	int i;
	int *m = (int*)malloc(sizeof(int)*l*c);
	int **M = (int**)malloc(sizeof(int*)*l);
	
	for(i=0;i<l;i++)
		M[i] = &m[c*i];
	
	return M;		
}

int main(int argc, char *argv[]){

	int i,j,k,v1,v2,p;	// auxiliares de loop e leitura de dados

	// _MPI_INIT_
	MPI_Status status;
	MPI_Init (&argc, &argv);
	MPI_Comm_rank(MPI_COMM_WORLD, &myrank);	
	MPI_Comm_size(MPI_COMM_WORLD, &P); 

	// Numero de threads a serem criadas em cada processo
	omp_set_num_threads(NUM_THREADS);

	// Le dimensoes das matrizes de entrada
	if(myrank==0)scanf("%d%d",&V,&E);

	// Envia dimensoes para nos escravos
	MPI_Bcast(&V, 1, MPI_INT, 0, MPI_COMM_WORLD);
	MPI_Bcast(&E, 1, MPI_INT, 0, MPI_COMM_WORLD);

	// Aloca matriz de distâncias e uma auxiliar.
	D = allocMatrix(V,V);
	// Aloca lista de arestas
	edge = allocMatrix(E,3);

	// No mestre le os valores de entrada das matrizes A e B
	if(myrank==0){
		// Inicia matriz de distancias
		for(i=0;i<V;i++){
			for(j=0;j<V;j++){
				if(i==j)
					D[i][j] = 0;
				else
					D[i][j] = INF;
			}
		}

		// Le matriz de pesos e lista de arestas
		for(i=0;i<E;i++){
			scanf("%d%d%d",&v1,&v2,&p);
			
			D[v1-1][v2-1] = p;
			D[v2-1][v1-1] = p;

			edge[i][0] = v1-1;
			edge[i][1] = v2-1;
			edge[i][2] = p;
		}

	}

	// Envia matriz D inteira para os nos escravos
	MPI_Bcast(D[0],V*V,MPI_INT, 0, MPI_COMM_WORLD);
	// Envia lista de arestas para os nos escravos
	MPI_Bcast(edge[0],3*E,MPI_INT, 0, MPI_COMM_WORLD);

	// Calcula a matriz de distancias minimas do grafo
	runFloydWarshall();

	// Procura o centro absoluto do grafo e retorna o valor de seu raio.
	findAbsoluteCenter();

	if(myrank==0) printf("%.2f\n",L);

	free(D);
	free(edge);
	
	// Fim.
	MPI_Finalize();
	
	return 0;
}

